! modal_aero_gasaerexch.F90


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!BOP
!
! !MODULE: modal_aero_gasaerexch --- does modal aerosol gas-aerosol exchange
!
! !INTERFACE:
   module modal_aero_gasaerexch

! !USES:
  use shr_kind_mod,    only:  r8 => shr_kind_r8
  use chem_mods,       only:  gas_pcnst
  use modal_aero_data, only:  maxd_aspectype

  implicit none
  private
  save

! !PUBLIC MEMBER FUNCTIONS:
  public modal_aero_gasaerexch_sub, modal_aero_gasaerexch_init

! !PUBLIC DATA MEMBERS:
  integer, parameter :: pcnstxx = gas_pcnst
  integer, parameter, public :: maxspec_pcage = maxd_aspectype

  integer, public :: npair_pcage
  integer, public :: modefrm_pcage
  integer, public :: modetoo_pcage
  integer, public :: nspecfrm_pcage
  integer, public :: lspecfrm_pcage(maxspec_pcage)
  integer, public :: lspectoo_pcage(maxspec_pcage)


! real(r8), parameter, public :: n_so4_monolayers_pcage = 1.0_r8
  real(r8), parameter, public :: n_so4_monolayers_pcage = 3.0_r8
! number of so4(+nh4) monolayers needed to "age" a carbon particle

  real(r8), parameter, public :: &
              dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10
! thickness of the so4 monolayers (m)
! for so4(+nh4), use bi-sulfate mw and 1.77 g/cm3,
!    --> 1 mol so4(+nh4)  = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
! aging criterion is approximate so do not try to distinguish
!    sulfuric acid, bisulfate, ammonium sulfate

  real(r8), save, public :: soa_equivso4_factor
! this factor converts an soa volume to a volume of so4(+nh4)
! having same hygroscopicity as the soa


! !DESCRIPTION: This module implements ...
!
! !REVISION HISTORY:
!
!   RCE 07.04.13:  Adapted from MIRAGE2 code
!
!EOP
!----------------------------------------------------------------------
!BOC

! list private module data here

!EOC
!----------------------------------------------------------------------


  contains


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!BOP
! !ROUTINE:  modal_aero_gasaerexch_sub --- ...
!
! !INTERFACE:
subroutine modal_aero_gasaerexch_sub(                            &
                        lchnk,    ncol,     nstep,               &
                        loffset,  deltat,                        &
                        latndx,   lonndx,                        &
                        t,        pmid,     pdel,                &
                        q,                  qqcw,                &
                        dqdt_other,         dqqcwdt_other,       &
                        dgncur_a,           dgncur_awet          )

! !USES:
use modal_aero_data
use modal_aero_rename, only:  modal_aero_rename_sub

use cam_history,       only:  outfld, fieldname_len
use chem_mods,         only:  adv_mass
use constituents,      only:  pcnst, cnst_name, cnst_get_ind
use mo_tracname,       only:  solsym
use physconst,         only:  gravit, mwdry, rair
use ppgrid,            only:  pcols, pver
                                                                                                                                            
use abortutils,        only : endrun
use spmd_utils,        only : iam, masterproc
                                                                                                                                            


implicit none

! !PARAMETERS:
   integer,  intent(in)    :: lchnk                ! chunk identifier
   integer,  intent(in)    :: ncol                 ! number of atmospheric column
   integer,  intent(in)    :: nstep                ! model time-step number
   integer,  intent(in)    :: loffset              ! offset applied to modal aero "ptrs"
   integer,  intent(in)    :: latndx(pcols), lonndx(pcols)
   real(r8), intent(in)    :: deltat               ! time step (s)

   real(r8), intent(inout) :: q(ncol,pver,pcnstxx) ! tracer mixing ratio (TMR) array
                                                   ! *** MUST BE  #/kmol-air for number
                                                   ! *** MUST BE mol/mol-air for mass
                                                   ! *** NOTE ncol dimension
   real(r8), intent(inout) :: qqcw(ncol,pver,pcnstxx) 
                                                   ! like q but for cloud-borner tracers
   real(r8), intent(in)    :: dqdt_other(ncol,pver,pcnstxx) 
                                                   ! TMR tendency from other continuous
                                                   ! growth processes (aqchem, soa??)
                                                   ! *** NOTE ncol dimension
   real(r8), intent(in)    :: dqqcwdt_other(ncol,pver,pcnstxx) 
                                                   ! like dqdt_other but for cloud-borner tracers
   real(r8), intent(in)    :: t(pcols,pver)        ! temperature at model levels (K)
   real(r8), intent(in)    :: pmid(pcols,pver)     ! pressure at model levels (Pa)
   real(r8), intent(in)    :: pdel(pcols,pver)     ! pressure thickness of levels (Pa)
   real(r8), intent(in)    :: dgncur_a(pcols,pver,ntot_amode)
   real(r8), intent(in)    :: dgncur_awet(pcols,pver,ntot_amode)
                                 ! dry & wet geo. mean dia. (m) of number distrib.

! !DESCRIPTION: 
! computes TMR (tracer mixing ratio) tendencies for gas condensation
!    onto aerosol particles
!
! this version does condensation of H2SO4, NH3, and MSA, both treated as
! completely non-volatile (gas --> aerosol, but no aerosol --> gas)
!    gas H2SO4 goes to aerosol SO4
!    gas MSA (if present) goes to aerosol SO4
!       aerosol MSA is not distinguished from aerosol SO4
!    gas NH3 (if present) goes to aerosol NH4
!       if gas NH3 is not present, then ????
!  
!
! !REVISION HISTORY:
!   RCE 07.04.13:  Adapted from MIRAGE2 code
!
!EOP
!----------------------------------------------------------------------
!BOC

! local variables
   integer, parameter :: jsrflx_gaexch = 1
   integer, parameter :: jsrflx_rename = 2
   integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
   integer, parameter :: method_soa = 2
!     method_soa=0 is no uptake
!     method_soa=1 is irreversible uptake done like h2so4 uptake
!     method_soa=2 is reversible uptake using subr modal_aero_soaexch

   integer :: i, icol_diag, iq
   integer :: idiagss
   integer :: ido_so4a(ntot_amode), ido_nh4a(ntot_amode), ido_soaa(ntot_amode)
   integer :: jac, jsrf
   integer :: k
   integer :: l, l2, lb, lsfrm, lstoo
   integer :: l_so4g, l_nh4g, l_msag, l_soag
   integer :: n, niter, niter_max, ntot_soamode

   logical :: is_dorename_atik, dorename_atik(ncol,pver)

   character(len=fieldname_len+3) :: fieldname

   real (r8) :: avg_uprt_nh4, avg_uprt_so4, avg_uprt_soa
   real (r8) :: deltatxx
   real (r8) :: dqdt_nh4(ntot_amode), dqdt_so4(ntot_amode), dqdt_soa(ntot_amode)
   real (r8) :: fac_m2v_nh4, fac_m2v_so4, fac_m2v_soa
   real (r8) :: fac_m2v_pcarbon(maxd_aspectype)
   real (r8) :: fac_volsfc_pcarbon
   real (r8) :: fgain_nh4(ntot_amode), fgain_so4(ntot_amode), fgain_soa(ntot_amode)
   real (r8) :: g0_soa
   real (r8) :: pdel_fac
   real (r8) :: qmax_nh4, qnew_nh4, qnew_so4
   real (r8) :: qold_nh4(ntot_amode), qold_so4(ntot_amode)
   real (r8) :: qold_soa(ntot_amode), qold_poa(ntot_amode)
   real (r8) :: sum_dqdt_msa, sum_dqdt_so4, sum_dqdt_soa
   real (r8) :: sum_dqdt_nh4, sum_dqdt_nh4_b
   real (r8) :: sum_uprt_msa, sum_uprt_nh4, sum_uprt_so4, sum_uprt_soa
   real (r8) :: tmp1, tmp2
   real (r8) :: uptkrate(ntot_amode,pcols,pver)  
   real (r8) :: uptkratebb(ntot_amode), uptkrate_soa(ntot_amode)  
                ! gas-to-aerosol mass transfer rates (1/s)
   real (r8) :: vol_core, vol_shell
   real (r8) :: xferfrac_pcage, xferfrac_max
   real (r8) :: xferrate

   logical  :: do_msag         ! true if msa gas is a species
   logical  :: do_nh4g         ! true if nh3 gas is a species
   logical  :: do_soag         ! true if soa gas is a species

   logical  :: dotend(pcnstxx)          ! identifies species directly involved in
                                        !    gas-aerosol exchange (gas condensation)
   logical  :: dotendqqcw(pcnstxx)      ! like dotend but for cloud-borner tracers
   logical  :: dotendrn(pcnstxx), dotendqqcwrn(pcnstxx)
                                        ! identifies species involved in renaming
                                        !    after "continuous growth"
                                        !    (gas-aerosol exchange and aqchem)

   integer, parameter :: nsrflx = 2     ! last dimension of qsrflx
   real(r8) :: dqdt(ncol,pver,pcnstxx)  ! TMR "delta q" array - NOTE dims
   real(r8) :: dqqcwdt(ncol,pver,pcnstxx) ! like dqdt but for cloud-borner tracers
   real(r8) :: qsrflx(pcols,pcnstxx,nsrflx)
                              ! process-specific column tracer tendencies 
                              ! (1=renaming, 2=gas condensation)
   real(r8) :: qqcwsrflx(pcols,pcnstxx,nsrflx)

!  following only needed for diagnostics
   real(r8) :: qold(ncol,pver,pcnstxx)  ! NOTE dims
   real(r8) :: qnew(ncol,pver,pcnstxx)  ! NOTE dims
   real(r8) :: qdel(ncol,pver,pcnstxx)  ! NOTE dims
   real(r8) :: dumavec(1000), dumbvec(1000), dumcvec(1000)
   real(r8) :: qqcwold(ncol,pver,pcnstxx)
   real(r8) :: dqdtsv1(ncol,pver,pcnstxx)
   real(r8) :: dqqcwdtsv1(ncol,pver,pcnstxx)


!----------------------------------------------------------------------


   icol_diag = -1
   if (ldiag1 > 0) then
   if (nstep < 3) then
      do i = 1, ncol
         if ((latndx(i) == 23) .and. (lonndx(i) == 37)) icol_diag = i
      end do
   end if
   end if


! set gas species indices
   call cnst_get_ind( 'H2SO4', l_so4g, .false. )
   call cnst_get_ind( 'NH3',   l_nh4g, .false. )
   call cnst_get_ind( 'MSA',   l_msag, .false. )
   call cnst_get_ind( 'SOAG',  l_soag, .false. )
   l_so4g = l_so4g - loffset
   l_nh4g = l_nh4g - loffset
   l_msag = l_msag - loffset
   l_soag = l_soag - loffset
   if ((l_so4g <= 0) .or. (l_so4g > pcnstxx)) then
      write( *, '(/a/a,2i7)' )   &
         '*** modal_aero_gasaerexch_sub -- cannot find H2SO4 species',   &
         '    l_so4g, loffset =', l_so4g, loffset
      call endrun( 'modal_aero_gasaerexch_sub error' )
   end if
   do_nh4g = .false.
   do_msag = .false.
   if ((l_nh4g > 0) .and. (l_nh4g <= pcnstxx)) do_nh4g = .true.
   if ((l_msag > 0) .and. (l_msag <= pcnstxx)) do_msag = .true.
   do_soag = .false.
   if ((method_soa == 1) .or. (method_soa == 2)) then
      if ((l_soag > 0) .and. (l_soag <= pcnstxx)) do_soag = .true.
   else if (method_soa /= 0) then
      write(*,'(/a,1x,i10)') '*** modal_aero_gasaerexch_sub - bad method_soa =', method_soa
      call endrun( 'modal_aero_gasaerexch_sub error' )
   end if

! set tendency flags
   dotend(:) = .false.
   dotendqqcw(:) = .false.
   ido_so4a(:) = 0
   ido_nh4a(:) = 0
   ido_soaa(:) = 0

   dotend(l_so4g) = .true.
   if ( do_nh4g ) dotend(l_nh4g) = .true.
   if ( do_msag ) dotend(l_msag) = .true.
   if ( do_soag ) dotend(l_soag) = .true.
   ntot_soamode = 0
   do n = 1, ntot_amode
      l = lptr_so4_a_amode(n)-loffset
      if ((l > 0) .and. (l <= pcnstxx)) then
         dotend(l) = .true.
         ido_so4a(n) = 1
         if ( do_nh4g ) then
            l = lptr_nh4_a_amode(n)-loffset
            if ((l > 0) .and. (l <= pcnstxx)) then
               dotend(l) = .true.
               ido_nh4a(n) = 1
            end if
         end if
      end if
      if ( do_soag ) then
         l = lptr_soa_a_amode(n)-loffset
         if ((l > 0) .and. (l <= pcnstxx)) then
            dotend(l) = .true.
            ido_soaa(n) = 1
            ntot_soamode = n
         end if
      end if
   end do
   if ( do_soag ) ntot_soamode = max( ntot_soamode, modefrm_pcage )

   if (modefrm_pcage > 0) then
      ido_so4a(modefrm_pcage) = 2
      if (ido_nh4a(modetoo_pcage) == 1) ido_nh4a(modefrm_pcage) = 2
      if (ido_soaa(modetoo_pcage) == 1) ido_soaa(modefrm_pcage) = 2
      do iq = 1, nspecfrm_pcage
         lsfrm = lspecfrm_pcage(iq)-loffset
         lstoo = lspectoo_pcage(iq)-loffset
         if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
            dotend(lsfrm) = .true.
            if ((lstoo > 0) .and. (lstoo <= pcnst)) then
               dotend(lstoo) = .true.
            end if
         end if
      end do


      fac_m2v_so4 = specmw_so4_amode / specdens_so4_amode
      fac_m2v_nh4 = specmw_nh4_amode / specdens_nh4_amode
      fac_m2v_soa = specmw_soa_amode / specdens_soa_amode
      fac_m2v_pcarbon(:) = 0.0
      n = modeptr_pcarbon
      do l = 1, nspec_amode(n)
         l2 = lspectype_amode(l,n)
!   fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
!        [m3-AP/kmol-AP]    = [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
         fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
      end do
      fac_volsfc_pcarbon = exp( 2.5*(alnsg_amode(n)**2) )
      xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8)   ! 1-eps
   end if


! zero out tendencies and other
   dqdt(:,:,:) = 0.0
   dqqcwdt(:,:,:) = 0.0
   qsrflx(:,:,:) = 0.0
   qqcwsrflx(:,:,:) = 0.0

! compute gas-to-aerosol mass transfer rates
   call gas_aer_uptkrates( ncol,       loffset,                &
                           q,          t,          pmid,       &
                           dgncur_awet,            uptkrate    )


! use this for tendency calcs to avoid generating very small negative values
   deltatxx = deltat * (1.0d0 + 1.0d-15)


   jsrf = jsrflx_gaexch
   do k=1,pver
     do i=1,ncol

!   fgain_so4(n) = fraction of total h2so4 uptake going to mode n
!   fgain_nh4(n) = fraction of total  nh3  uptake going to mode n
        sum_uprt_so4 = 0.0
        sum_uprt_nh4 = 0.0
        sum_uprt_soa = 0.0
        do n = 1, ntot_amode
            uptkratebb(n) = uptkrate(n,i,k)
            if (ido_so4a(n) > 0) then
                fgain_so4(n) = uptkratebb(n)
                sum_uprt_so4 = sum_uprt_so4 + fgain_so4(n)
                if (ido_so4a(n) == 1) then
                    qold_so4(n) = q(i,k,lptr_so4_a_amode(n)-loffset)
                else
                    qold_so4(n) = 0.0
                end if
            else
                fgain_so4(n) = 0.0
                qold_so4(n) = 0.0
            end if

            if (ido_nh4a(n) > 0) then
!   2.08 factor is for gas diffusivity (nh3/h2so4)
!   differences in fuch-sutugin and accom coef ignored
                fgain_nh4(n) = uptkratebb(n)*2.08
                sum_uprt_nh4 = sum_uprt_nh4 + fgain_nh4(n)
                if (ido_nh4a(n) == 1) then
                    qold_nh4(n) = q(i,k,lptr_nh4_a_amode(n)-loffset)
                else
                    qold_nh4(n) = 0.0
                end if
            else
                fgain_nh4(n) = 0.0
                qold_nh4(n) = 0.0
            end if

            if (ido_soaa(n) > 0) then
!   0.81 factor is for gas diffusivity (soa/h2so4)
!   (differences in fuch-sutugin and accom coef ignored)
                fgain_soa(n) = uptkratebb(n)*0.81
                sum_uprt_soa = sum_uprt_soa + fgain_soa(n)
                if (ido_soaa(n) == 1) then
                    qold_soa(n) = q(i,k,lptr_soa_a_amode(n)-loffset)
                    l = lptr_pom_a_amode(n)-loffset
                    if (l > 0) then
                       qold_poa(n) = q(i,k,l)
                    else
                       qold_poa(n) = 0.0
                    end if
                else
                    qold_soa(n) = 0.0
                    qold_poa(n) = 0.0
                end if
            else
                fgain_soa(n) = 0.0
                qold_soa(n) = 0.0
                qold_poa(n) = 0.0
            end if
            uptkrate_soa(n) = fgain_soa(n)
        end do

        if (sum_uprt_so4 > 0.0) then
            do n = 1, ntot_amode
                fgain_so4(n) = fgain_so4(n) / sum_uprt_so4
            end do
        end if
!       at this point (sum_uprt_so4 <= 0.0) only when all the fgain_so4 are zero
        if (sum_uprt_nh4 > 0.0) then
            do n = 1, ntot_amode
                fgain_nh4(n) = fgain_nh4(n) / sum_uprt_nh4
            end do
        end if
        if (sum_uprt_soa > 0.0) then
            do n = 1, ntot_amode
                fgain_soa(n) = fgain_soa(n) / sum_uprt_soa
            end do
        end if

!   uptake amount (fraction of gas uptaken) over deltat
        avg_uprt_so4 = (1.0 - exp(-deltatxx*sum_uprt_so4))/deltatxx
        avg_uprt_nh4 = (1.0 - exp(-deltatxx*sum_uprt_nh4))/deltatxx
        avg_uprt_soa = (1.0 - exp(-deltatxx*sum_uprt_soa))/deltatxx
 
!   sum_dqdt_so4 = so4_a tendency from h2so4 gas uptake (mol/mol/s)
!   sum_dqdt_msa = msa_a tendency from msa   gas uptake (mol/mol/s)
!   sum_dqdt_nh4 = nh4_a tendency from nh3   gas uptake (mol/mol/s)
!   sum_dqdt_soa = soa_a tendency from soa   gas uptake (mol/mol/s)
        sum_dqdt_so4 = q(i,k,l_so4g) * avg_uprt_so4
        if ( do_msag ) then
            sum_dqdt_msa = q(i,k,l_msag) * avg_uprt_so4
        else
            sum_dqdt_msa = 0.0
        end if
        if ( do_nh4g ) then
            sum_dqdt_nh4 = q(i,k,l_nh4g) * avg_uprt_nh4
        else
            sum_dqdt_nh4 = 0.0
        end if
        if ( do_soag ) then
            sum_dqdt_soa = q(i,k,l_soag) * avg_uprt_soa
        else
            sum_dqdt_soa = 0.0
        end if

!   compute TMR tendencies for so4, nh4, msa interstial aerosol
!   due to simple gas uptake
        pdel_fac = pdel(i,k)/gravit
        sum_dqdt_nh4_b = 0.0
        do n = 1, ntot_amode
            dqdt_so4(n) = fgain_so4(n)*(sum_dqdt_so4 + sum_dqdt_msa)
 
            if ( do_nh4g ) then
                dqdt_nh4(n) = fgain_nh4(n)*sum_dqdt_nh4
                qnew_nh4 = qold_nh4(n) + dqdt_nh4(n)*deltat
                qnew_so4 = qold_so4(n) + dqdt_so4(n)*deltat
                qmax_nh4 = 2.0*qnew_so4
                if (qnew_nh4 > qmax_nh4) then
                    dqdt_nh4(n) = (qmax_nh4 - qold_nh4(n))/deltatxx
                end if
                sum_dqdt_nh4_b = sum_dqdt_nh4_b + dqdt_nh4(n)
            end if
        end do

        if (( do_soag ) .and. (method_soa > 1)) then
!   compute TMR tendencies for soag and soa interstial aerosol
!   using soa parameterization
           niter_max = 1000
           dqdt_soa(:) = 0.0

!          idiagss = 0
!          if (ldiag2 > 0) then
!          if ((i == icol_diag) .and. (mod(k-1,5) == 0)) then
!             idiagss = 1
!          end if
!          end if

!          call modal_aero_soaexch( dtfull, temp, pres, &
!             niter, niter_max, ntot_soamode, &
!             g_soa_in, a_soa_in, a_poa_in, xferrate, &
!             g_soa_tend, a_soa_tend )
           call modal_aero_soaexch( deltat, t(i,k), pmid(i,k), &
              niter, niter_max, ntot_soamode, &
              q(i,k,l_soag), qold_soa, qold_poa, uptkrate_soa, &
              tmp1, dqdt_soa )
!             tmp1, dqdt_soa, g0_soa, idiagss )
           sum_dqdt_soa = -tmp1

        else if ( do_soag ) then
!   compute TMR tendencies for soa interstial aerosol
!   due to simple gas uptake
           do n = 1, ntot_amode
                dqdt_soa(n) = fgain_soa(n)*sum_dqdt_soa
           end do
        else
           dqdt_soa(:) = 0.0
        end if
 
        do n = 1, ntot_amode
            if (ido_so4a(n) == 1) then
                l = lptr_so4_a_amode(n)-loffset
                dqdt(i,k,l) = dqdt_so4(n)
                qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(n)*pdel_fac
            end if
 
            if ( do_nh4g ) then
                if (ido_nh4a(n) == 1) then
                    l = lptr_nh4_a_amode(n)-loffset
                    dqdt(i,k,l) = dqdt_nh4(n)
                    qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(n)*pdel_fac
                end if
            end if
 
            if ( do_soag ) then
                if (ido_soaa(n) == 1) then
                    l = lptr_soa_a_amode(n)-loffset
                    dqdt(i,k,l) = dqdt_soa(n)
                    qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(n)*pdel_fac
                end if
            end if
        end do
 
!   compute TMR tendencies for h2so4, nh3, and msa gas
!   due to simple gas uptake
        l = l_so4g
        dqdt(i,k,l) = -sum_dqdt_so4
        qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac

        if ( do_msag ) then
           l = l_msag
           dqdt(i,k,l) = -sum_dqdt_msa
           qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
        end if

        if ( do_nh4g ) then
           l = l_nh4g
           dqdt(i,k,l) = -sum_dqdt_nh4_b
           qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
        end if
 
        if ( do_soag ) then
           l = l_soag
           dqdt(i,k,l) = -sum_dqdt_soa
           qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt(i,k,l)*pdel_fac
        end if
 
!   compute TMR tendencies associated with primary carbon aging
        if (modefrm_pcage > 0) then
           n = modeptr_pcarbon
           vol_shell = deltat *   &
                     ( dqdt_so4(n)*fac_m2v_so4 + dqdt_nh4(n)*fac_m2v_nh4 +   &
                       dqdt_soa(n)*fac_m2v_soa*soa_equivso4_factor )
           vol_core = 0.0
           do l = 1, nspec_amode(n)
              vol_core = vol_core + &
                    q(i,k,lmassptr_amode(l,n)-loffset)*fac_m2v_pcarbon(l)
           end do
!   ratio1 = vol_shell/vol_core = 
!      actual hygroscopic-shell-volume/carbon-core-volume after gas uptake
!   ratio2 = 6.0_r8*dr_so4_monolayers_pcage/(dgncur_a*fac_volsfc_pcarbon)
!      = (shell-volume corresponding to n_so4_monolayers_pcage)/core-volume 
!      The 6.0/(dgncur_a*fac_volsfc_pcarbon) = (mode-surface-area/mode-volume)
!   Note that vol_shell includes both so4+nh4 AND soa as "equivalent so4",
!      The soa_equivso4_factor accounts for the lower hygroscopicity of soa.
!
!   Define xferfrac_pcage = min( 1.0, ratio1/ratio2)
!   But ratio1/ratio2 == tmp1/tmp2, and coding below avoids possible overflow 
!
           tmp1 = vol_shell*dgncur_a(i,k,n)*fac_volsfc_pcarbon
           tmp2 = max( 6.0_r8*dr_so4_monolayers_pcage*vol_core, 0.0_r8 )
           if (tmp1 >= tmp2) then
              xferfrac_pcage = xferfrac_max
           else
              xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
           end if

           if (xferfrac_pcage > 0.0_r8) then
              do iq = 1, nspecfrm_pcage
                 lsfrm = lspecfrm_pcage(iq)-loffset
                 lstoo = lspectoo_pcage(iq)-loffset
                 xferrate = (xferfrac_pcage/deltat)*q(i,k,lsfrm)
                 dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferrate
                 qsrflx(i,lsfrm,jsrf) = qsrflx(i,lsfrm,jsrf) - xferrate*pdel_fac
                 if ((lstoo > 0) .and. (lstoo <= pcnst)) then
                     dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferrate
                     qsrflx(i,lstoo,jsrf) = qsrflx(i,lstoo,jsrf) + xferrate*pdel_fac
                 end if
              end do

              if (ido_so4a(modetoo_pcage) > 0) then
                 l = lptr_so4_a_amode(modetoo_pcage)-loffset
                 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_so4(modefrm_pcage)
                 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_so4(modefrm_pcage)*pdel_fac
              end if

              if (ido_nh4a(modetoo_pcage) > 0) then
                 l = lptr_nh4_a_amode(modetoo_pcage)-loffset
                 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_nh4(modefrm_pcage)
                 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_nh4(modefrm_pcage)*pdel_fac
              end if

              if (ido_soaa(modetoo_pcage) > 0) then
                 l = lptr_soa_a_amode(modetoo_pcage)-loffset
                 dqdt(i,k,l) = dqdt(i,k,l) + dqdt_soa(modefrm_pcage)
                 qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf) + dqdt_soa(modefrm_pcage)*pdel_fac
              end if
           end if

        end if


! diagnostics start -------------------------------------------------------
   if (ldiag2 > 0) then
   if (i == icol_diag) then 
   if (mod(k-1,5) == 0) then
      write(*,'(a,43i5)') 'gasaerexch aaa nstep,lat,lon,k', nstep, latndx(i), lonndx(i), k
      write(*,'(a,1p,10e12.4)') 'uptkratebb   ', uptkratebb(:)
      write(*,'(a,1p,10e12.4)') 'sum_uprt_so4 ', sum_uprt_so4
      write(*,'(a,1p,10e12.4)') 'fgain_so4    ', fgain_so4(:)
      write(*,'(a,1p,10e12.4)') 'sum_uprt_nh4 ', sum_uprt_nh4
      write(*,'(a,1p,10e12.4)') 'fgain_nh4    ', fgain_nh4(:)
      write(*,'(a,1p,10e12.4)') 'sum_uprt_soa ', sum_uprt_soa
      write(*,'(a,1p,10e12.4)') 'fgain_soa    ', fgain_soa(:)
      write(*,'(a,1p,10e12.4)') 'so4g o,dqdt,n', q(i,k,l_so4g), sum_dqdt_so4, &
                                                (q(i,k,l_so4g)-deltat*sum_dqdt_so4)
      write(*,'(a,1p,10e12.4)') 'nh3g o,dqdt,n', q(i,k,l_nh4g), sum_dqdt_nh4, sum_dqdt_nh4_b, &
                                                (q(i,k,l_nh4g)-deltat*sum_dqdt_nh4_b)
      write(*,'(a,1p,10e12.4)') 'soag o,dqdt,n', q(i,k,l_soag), sum_dqdt_soa, &
                                                (q(i,k,l_soag)-deltat*sum_dqdt_soa)
      write(*,'(a,i12,1p,10e12.4)') &
                                'method,g0,t,p', method_soa, g0_soa, t(i,k), pmid(i,k)
      write(*,'(a,1p,10e12.4)') 'so4 old      ', qold_so4(:)
      write(*,'(a,1p,10e12.4)') 'so4 dqdt     ', dqdt_so4(:)
      write(*,'(a,1p,10e12.4)') 'so4 new      ', (qold_so4(:)+deltat*dqdt_so4(:))
      write(*,'(a,1p,10e12.4)') 'nh4 old      ', qold_nh4(:)
      write(*,'(a,1p,10e12.4)') 'nh4 dqdt     ', dqdt_nh4(:)
      write(*,'(a,1p,10e12.4)') 'nh4 new      ', (qold_nh4(:)+deltat*dqdt_nh4(:))
      write(*,'(a,1p,10e12.4)') 'soa old      ', qold_soa(:)
      write(*,'(a,1p,10e12.4)') 'soa dqdt     ', dqdt_soa(:)
      write(*,'(a,1p,10e12.4)') 'soa new      ', (qold_soa(:)+deltat*dqdt_soa(:))
      write(*,'(a,1p,10e12.4)') 'vshell, core ', vol_shell, vol_core
      write(*,'(a,1p,10e12.4)') 'dr_mono, ... ', dr_so4_monolayers_pcage,   &
                                 soa_equivso4_factor
      write(*,'(a,1p,10e12.4)') 'dgn, ...     ', dgncur_a(i,k,modefrm_pcage),   &
                                 fac_volsfc_pcarbon
      write(*,'(a,1p,10e12.4)') 'tmp1, tmp2   ', tmp1, tmp2
      write(*,'(a,1p,10e12.4)') 'xferfrac_age ', xferfrac_pcage
   end if
   end if
   end if
! diagnostics end ---------------------------------------------------------


     end do   ! "i = 1, ncol"
   end do     ! "k = 1, pver"


! set "temporary testing arrays"
   qold(:,:,:) = q(:,:,:)
   qqcwold(:,:,:) = qqcw(:,:,:)
   dqdtsv1(:,:,:) = dqdt(:,:,:)
   dqqcwdtsv1(:,:,:) = dqqcwdt(:,:,:)


!
! do renaming calcs
!
   dotendrn(:) = .false.
   dotendqqcwrn(:) = .false.
   dorename_atik(1:ncol,:) = .true.
   is_dorename_atik = .true.
   if (ncol >= -13579) then
      call modal_aero_rename_sub(                              &
                       'modal_aero_gasaerexch_sub',            &
                       lchnk,             ncol,      nstep,    &
                       loffset,           deltat,              &
                       latndx,            lonndx,              &
                       pdel,                                   &
                       dotendrn,          q,                   &
                       dqdt,              dqdt_other,          &
                       dotendqqcwrn,      qqcw,                &
                       dqqcwdt,           dqqcwdt_other,       &
                       is_dorename_atik,  dorename_atik,       &
                       jsrflx_rename,     nsrflx,              &
                       qsrflx,            qqcwsrflx            )
   end if


!
!  apply the dqdt to update q (and same for qqcw)
!
   do l = 1, pcnstxx
      if ( dotend(l) .or. dotendrn(l) ) then
         do k = 1, pver
         do i = 1, ncol
            q(i,k,l) = q(i,k,l) + dqdt(i,k,l)*deltat
         end do
         end do
      end if
      if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
         do k = 1, pver
         do i = 1, ncol
            qqcw(i,k,l) = qqcw(i,k,l) + dqqcwdt(i,k,l)*deltat
         end do
         end do
      end if
   end do


! diagnostics start -------------------------------------------------------
   if (ldiag3 > 0) then
   if (icol_diag > 0) then 
      i = icol_diag
      write(*,'(a,3i5)') 'gasaerexch ppp nstep,lat,lon', nstep, latndx(i), lonndx(i)
      write(*,'(2i5,3(2x,a))') 0, 0, 'ppp', 'pdel for all k'
      write(*,'(1p,7e12.4)') (pdel(i,k), k=1,pver)

      write(*,'(a,3i5)') 'gasaerexch ddd nstep,lat,lon', nstep, latndx(i), lonndx(i)
      do l = 1, pcnstxx
         lb = l + loffset

         if ( dotend(l) .or. dotendrn(l) ) then
            write(*,'(2i5,3(2x,a))') 1, l, 'ddd1', cnst_name(lb),    'qold for all k'
            write(*,'(1p,7e12.4)') (qold(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 1, l, 'ddd2', cnst_name(lb),    'qnew for all k'
            write(*,'(1p,7e12.4)') (q(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 1, l, 'ddd3', cnst_name(lb),    'dqdt from conden for all k'
            write(*,'(1p,7e12.4)') (dqdtsv1(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 1, l, 'ddd4', cnst_name(lb),    'dqdt from rename for all k'
            write(*,'(1p,7e12.4)') ((dqdt(i,k,l)-dqdtsv1(i,k,l)), k=1,pver)
            write(*,'(2i5,3(2x,a))') 1, l, 'ddd5', cnst_name(lb),    'dqdt other for all k'
            write(*,'(1p,7e12.4)') (dqdt_other(i,k,l), k=1,pver)
         end if

         if ( dotendqqcw(l) .or. dotendqqcwrn(l) ) then
            write(*,'(2i5,3(2x,a))') 2, l, 'ddd1', cnst_name_cw(lb), 'qold for all k'
            write(*,'(1p,7e12.4)') (qqcwold(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 2, l, 'ddd2', cnst_name_cw(lb), 'qnew for all k'
            write(*,'(1p,7e12.4)') (qqcw(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 2, l, 'ddd3', cnst_name_cw(lb), 'dqdt from conden for all k'
            write(*,'(1p,7e12.4)') (dqqcwdtsv1(i,k,l), k=1,pver)
            write(*,'(2i5,3(2x,a))') 2, l, 'ddd4', cnst_name_cw(lb), 'dqdt from rename for all k'
            write(*,'(1p,7e12.4)') ((dqqcwdt(i,k,l)-dqqcwdtsv1(i,k,l)), k=1,pver)
            write(*,'(2i5,3(2x,a))') 2, l, 'ddd5', cnst_name_cw(lb), 'dqdt other for all k'
            write(*,'(1p,7e12.4)') (dqqcwdt_other(i,k,l), k=1,pver)
         end if

      end do

      write(*,'(a,3i5)') 'gasaerexch fff nstep,lat,lon', nstep, latndx(i), lonndx(i)
      do l = 1, pcnstxx
         lb = l + loffset
         if ( dotend(l) .or. dotendrn(l) .or. dotendqqcw(l) .or. dotendqqcwrn(l) ) then
            write(*,'(i5,2(2x,a,2l3))') l, &
               cnst_name(lb), dotend(l), dotendrn(l), &
               cnst_name_cw(lb), dotendqqcw(l), dotendqqcwrn(l)
         end if
      end do

   end if
   end if
! diagnostics end ---------------------------------------------------------


!   do history file column-tendency fields
	do l = 1, pcnstxx
	lb = l + loffset

	do jsrf = 1, 2

	do jac = 1, 2

	    if (jac == 1) then
		if (jsrf == jsrflx_gaexch) then
		    if ( .not. dotend(l) ) cycle
		    fieldname = trim(cnst_name(lb)) // '_sfgaex1'
		else if (jsrf == jsrflx_rename) then
		    if ( .not. dotendrn(l) ) cycle
		    fieldname = trim(cnst_name(lb)) // '_sfgaex2'
		else 
		    cycle
		end if
		do i = 1, ncol
		    qsrflx(i,l,jsrf) = qsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
		end do
		call outfld( fieldname, qsrflx(:,l,jsrf), pcols, lchnk )

	    else
		if (jsrf == jsrflx_gaexch) then
		    cycle
		else if (jsrf == jsrflx_rename) then
		    if ( .not. dotendqqcwrn(l) ) cycle
		    fieldname = trim(cnst_name_cw(lb)) // '_sfgaex2'
		else 
		    cycle
		end if
		do i = 1, ncol
		    qqcwsrflx(i,l,jsrf) = qqcwsrflx(i,l,jsrf)*(adv_mass(l)/mwdry)
		end do
		call outfld( fieldname, qqcwsrflx(:,l,jsrf), pcols, lchnk )
	    end if

!	    if (( masterproc ) .and. (nstep < 1)) &
!		write(*,'(2(a,2x),1p,e11.3)') &
!		'modal_aero_newnuc_sub outfld', fieldname, adv_mass(l)

 	    if (ldiag4 > 0) then
 	    if (icol_diag > 0) then
 		i = icol_diag
 		if (jac == 1) then
 		    tmp1 = qsrflx(i,l,jsrf)
 		else
 		    tmp1 = qqcwsrflx(i,l,jsrf)
 		end if
 		write(*,'(a,4i5,2x,a,1p,2e12.4)')   &
 		    'gasaerexch nstep,lat,lon,l,fieldname,qsrflx,adv_mass',   &
 		    nstep, latndx(i), lonndx(i), l, fieldname, tmp1, adv_mass(l)
 	    end if
 	    end if

	end do ! jac = ...
	end do ! jsrf = ...
	end do ! l = ...


	return
!EOC
   end subroutine modal_aero_gasaerexch_sub


!----------------------------------------------------------------------
!----------------------------------------------------------------------
subroutine gas_aer_uptkrates( ncol,       loffset,                &
                              q,          t,          pmid,       &
                              dgncur_awet,            uptkrate    )

!
!                         /
!   computes   uptkrate = | dx  dN/dx  gas_conden_rate(Dp(x))
!                         /
!   using Gauss-Hermite quadrature of order nghq=2
!
!       Dp = particle diameter (cm)
!       x = ln(Dp)
!       dN/dx = log-normal particle number density distribution
!       gas_conden_rate(Dp) = 2 * pi * gasdiffus * Dp * F(Kn,ac)
!           F(Kn,ac) = Fuchs-Sutugin correction factor
!           Kn = Knudsen number
!           ac = accomodation coefficient
!

!use modal_aero_data, only:  ntot_amode, ntot_amode, nspec_amode,   &
!                            lspectype_amode, lmassptr_amode,   &
!                            sigmag_amode,   &
!                            specdens_amode, specmw_amode
use modal_aero_data, only:  ntot_amode, ntot_amode,   &
                            numptr_amode,   &
                            sigmag_amode

use ppgrid
use constituents, only: pcnst, cnst_name
use physconst, only: mwdry, rair

implicit none


   integer,  intent(in) :: ncol                 ! number of atmospheric column
   integer,  intent(in) :: loffset
   real(r8), intent(in) :: q(ncol,pver,pcnstxx) ! Tracer array (mol,#/mol-air)
   real(r8), intent(in) :: t(pcols,pver)        ! Temperature in Kelvin
   real(r8), intent(in) :: pmid(pcols,pver)     ! Air pressure in Pa
   real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)

   real(r8), intent(out) :: uptkrate(ntot_amode,pcols,pver)  
                            ! gas-to-aerosol mass transfer rates (1/s)


! local 
   integer, parameter :: nghq = 2
   integer :: i, iq, k, l1, l2, la, n

   real(r8), parameter :: tworootpi = 3.5449077
   real(r8), parameter :: root2 = 1.4142135
   real(r8), parameter :: beta = 2.0

   real(r8) :: aircon
   real(r8) :: const
   real(r8) :: dp, dum_m2v
   real(r8) :: dryvol_a(pcols,pver)
   real(r8) :: gasdiffus, gasspeed
   real(r8) :: freepathx2, fuchs_sutugin
   real(r8) :: knudsen
   real(r8) :: lndp, lndpgn, lnsg
   real(r8) :: num_a
   real(r8) :: rhoair
   real(r8) :: sumghq
   real(r8), save :: xghq(nghq), wghq(nghq) ! quadrature abscissae and weights

   data xghq / 0.70710678, -0.70710678 /
   data wghq / 0.88622693,  0.88622693 /


! outermost loop over all modes
   do n = 1, ntot_amode

! 22-aug-2007 rc easter - get number from q array rather
!    than computing a "bounded" number conc.
!! compute dry volume = sum_over_components{ component_mass / density }
!!    (m3-AP/mol-air)
!! compute it for all i,k to improve accessing q array
!      dryvol_a(1:ncol,:) = 0.0_r8
!      do l1 = 1, nspec_amode(n)
!         l2 = lspectype_amode(l1,n)
!! dum_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
!! [m3-AP/kmol-AP]= [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
!         dum_m2v = specmw_amode(l2) / specdens_amode(l2)
!         la = lmassptr_amode(l1,n)
!         dryvol_a(1:ncol,:) = dryvol_a(1:ncol,:)    &
!                            + max(0.0_r8,q(1:ncol,:,la))*dum_m2v
!      end do

! loops k and i
      do k=1,pver
      do i=1,ncol

         rhoair = pmid(i,k)/(rair*t(i,k))   ! (kg-air/m3)
!        aircon = 1.0e3*rhoair/mwdry        ! (mol-air/m3)

!!   "bounded" number conc. (#/m3)
!        num_a = dryvol_a(i,k)*v2ncur_a(i,k,n)*aircon

!   number conc. (#/m3) -- note q(i,k,numptr) is (#/kmol-air) 
!   so need aircon in (kmol-air/m3)
         aircon = rhoair/mwdry              ! (kmol-air/m3)
         num_a = q(i,k,numptr_amode(n)-loffset)*aircon

!   gasdiffus = h2so4 gas diffusivity from mosaic code (m^2/s)
!               (pmid must be Pa)
         gasdiffus = 0.557e-4 * (t(i,k)**1.75) / pmid(i,k)
!   gasspeed = h2so4 gas mean molecular speed from mosaic code (m/s)
         gasspeed  = 1.470e1 * sqrt(t(i,k))
!   freepathx2 = 2 * (h2so4 mean free path)  (m)
         freepathx2 = 6.0*gasdiffus/gasspeed

         lnsg   = log( sigmag_amode(n) )
         lndpgn = log( dgncur_awet(i,k,n) )   ! (m)
         const  = tworootpi * num_a * exp(beta*lndpgn + 0.5*(beta*lnsg)**2)
         
!   sum over gauss-hermite quadrature points
         sumghq = 0.0
         do iq = 1, nghq
            lndp = lndpgn + beta*lnsg**2 + root2*lnsg*xghq(iq)
            dp = exp(lndp)

!   knudsen number
            knudsen = freepathx2/dp
!   following assumes accomodation coefficient = ac = 0.65
!   (Adams & Seinfeld, 2002, JGR, and references therein)
!           fuchs_sutugin = (0.75*ac*(1. + knudsen)) /
!                           (knudsen*(1.0 + knudsen + 0.283*ac) + 0.75*ac)
            fuchs_sutugin = (0.4875*(1. + knudsen)) /   &
                            (knudsen*(1.184 + knudsen) + 0.4875)

            sumghq = sumghq + wghq(iq)*dp*fuchs_sutugin/(dp**beta)
         end do
         uptkrate(n,i,k) = const * gasdiffus * sumghq    

      end do   ! "do i = 1, ncol"
      end do   ! "do k = 1, pver"

   end do   ! "do n = 1, ntot_soamode"


   return
   end subroutine gas_aer_uptkrates

!----------------------------------------------------------------------


      subroutine modal_aero_soaexch( dtfull, temp, pres, &
          niter, niter_max, ntot_soamode, &
          g_soa_in, a_soa_in, a_poa_in, xferrate, &
          g_soa_tend, a_soa_tend )
!         g_soa_tend, a_soa_tend, g0_soa, idiagss )

!-----------------------------------------------------------------------
!
! Purpose:
!
! calculates condensation/evaporation of "soa gas" 
! to/from multiple aerosol modes in 1 grid cell
!
! key assumptions
! (1) ambient equilibrium vapor pressure of soa gas
!     is given by p0_soa_298 and delh_vap_soa
! (2) equilibrium vapor pressure of soa gas at aerosol
!     particle surface is given by raoults law in the form
!     g_star = g0_soa*[a_soa/(a_soa + a_opoa)]
! (3) (oxidized poa)/(total poa) is equal to frac_opoa (constant)
! 
!
! Author: R. Easter and R. Zaveri
!
!-----------------------------------------------------------------------
      implicit none

      real(r8), intent(in)  :: dtfull     ! full integration time step (s)
      real(r8), intent(in)  :: temp       ! air temperature (K)
      real(r8), intent(in)  :: pres       ! air pressure (Pa)
      integer,  intent(out) :: niter      ! number of iterations performed
      integer,  intent(in)  :: niter_max  ! max allowed number of iterations
      integer,  intent(in)  :: ntot_soamode             ! number of modes having soa
      real(r8), intent(in)  :: g_soa_in                 ! initial soa gas mixrat (mol/mol)
      real(r8), intent(in)  :: a_soa_in(ntot_soamode)   ! initial soa aerosol mixrat (mol/mol)
      real(r8), intent(in)  :: a_poa_in(ntot_soamode)   ! initial poa aerosol mixrat (mol/mol)
      real(r8), intent(in)  :: xferrate(ntot_soamode)   ! gas-aerosol mass transfer rate (1/s)
      real(r8), intent(out) :: g_soa_tend               ! soa gas mixrat tendency (mol/mol/s)
      real(r8), intent(out) :: a_soa_tend(ntot_soamode) ! soa aerosol mixrat tendency (mol/mol/s)
!     real(r8), intent(out) :: g0_soa   ! ambient soa gas equilib mixrat (mol/mol)
!     integer,  intent(in)  :: idiagss

      integer :: luna=6
      integer :: m

      real(r8), parameter :: alpha = 0.05  ! parameter used in calc of time step
      real(r8), parameter :: g_min1 = 1.0e-20
      real(r8), parameter :: opoa_frac = 0.1  ! fraction of poa that is opoa
      real(r8), parameter :: delh_vap_soa = 156.0e3
      ! delh_vap_soa = heat of vaporization for gas soa (J/mol)
      real(r8), parameter :: p0_soa_298 = 1.0e-10
      ! p0_soa_298 = soa gas equilib vapor presssure (atm) at 298 k
      real(r8), parameter :: rgas = 8.3144   ! gas constant in J/K/mol

      real(r8) :: a_opoa(ntot_soamode)    ! oxidized-poa aerosol mixrat (mol/mol)
      real(r8) :: a_soa(ntot_soamode)     ! soa aerosol mixrat (mol/mol)
      real(r8) :: a_soa_tmp               ! temporary soa aerosol mixrat (mol/mol)
      real(r8) :: beta(ntot_soamode)      ! dtcur*xferrate
      real(r8) :: dtcur    ! current time step (s)
      real(r8) :: dtmax    ! = (dtfull-tcur)
      real(r8) :: g_soa    ! soa gas mixrat (mol/mol)
      real(r8) :: g0_soa   ! ambient soa gas equilib mixrat (mol/mol)
      real(r8) :: g_star(ntot_soamode)    ! soa gas mixrat that is in equilib
                                          ! with each aerosol mode (mol/mol)
      real(r8) :: phi(ntot_soamode)       ! "relative driving force"
      real(r8) :: p0_soa   ! soa gas equilib vapor presssure (atm)
      real(r8) :: sat(ntot_soamode)
      real(r8) :: tcur     ! current integration time (from 0 s)
      real(r8) :: tmpa, tmpb
      real(r8) :: tot_soa  ! g_soa + sum( a_soa(:) )


! force things to be non-negative and calc tot_soa
! calc a_opoa (always slightly >0)
      g_soa = max( g_soa_in, 0.0_r8 )
      tot_soa = g_soa
      do m = 1, ntot_soamode
         a_soa(m) = max( a_soa_in(m), 0.0_r8 )
         tot_soa = tot_soa + a_soa(m)
         a_opoa(m) = opoa_frac*a_poa_in(m)
         a_opoa(m) = max( a_opoa(m), 1.0e-20_r8 )  ! force to small non-zero value
      end do

! calc ambient equilibrium soa gas
      p0_soa = p0_soa_298 * &
               exp( -(delh_vap_soa/rgas)*((1.0/temp)-(1.0/298.0)) )
      g0_soa = 1.01325e5*p0_soa/pres
! molecular weight adjustment
! the soa parameterization assumes that real gsoa and asoa have mw=150
! currently in cam3,
!    mw=12 is used (this has to do with the mozart preprocessor)
!    soag emission files (molec/cm2/s units) are set to give the desired 
!       mass emissions (kg/m2/s) and mass mixing ratios (kg/kg)
!       when mw=12 is applied
!    as a result, the molar mixing ratios for both gsoa and asoa
!       are artificially scaled up by (150/12)
!       and g0_soa must be similarly scaled up
      g0_soa = g0_soa*(150.0/12.0)
!     g0_soa = 0.0   ! force irreversible uptake

      niter = 0
      tcur = 0.0
      dtcur = 0.0
      phi(:) = 0.0
      g_star(:) = 0.0

!     if (idiagss > 0) then
!        write(luna,'(a,1p,10e11.3)') 'p0, g0_soa', p0_soa, g0_soa
!        write(luna,'(3a)') &
!           'niter, tcur,   dtcur,    phi(:),                       ', &
!           'g_star(:),                    ', &
!           'a_soa(:),                     g_soa'      
!        write(luna,'(3a)') &
!           '                         sat(:),                       ', &
!           'sat(:)*a_soa(:)               ', &
!           'a_opoa(:)'
!        write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
!           phi(:), g_star(:), a_soa(:), g_soa
!     end if


! integration loop -- does multiple substeps to reach dtfull
timeloop: do while (tcur < dtfull-1.0e-3_r8 )

      niter = niter + 1
      if (niter > niter_max) exit

      tmpa = 0.0
      do m = 1, ntot_soamode
         sat(m) = g0_soa/(a_soa(m) + a_opoa(m))
         g_star(m) = sat(m)*a_soa(m)
         phi(m) = (g_soa - g_star(m))/max(g_soa,g_star(m),g_min1)
         tmpa = tmpa + xferrate(m)*abs(phi(m))
      end do

      dtmax = dtfull-tcur
      if (dtmax*tmpa <= alpha) then
! here alpha/tmpa >= dtmax, so this is final substep
         dtcur = dtmax
         tcur = dtfull
      else
         dtcur = alpha/tmpa
         tcur = tcur + dtcur
      end if

! step 1 - for modes where soa is condensing, estimate "new" a_soa(m)
!    using an explicit calculation with "old" g_soa
!    and g_star(m) calculated using "old" a_soa(m)
! do this to get better estimate of "new" a_soa(m) and sat(m)
      do m = 1, ntot_soamode
         beta(m) = dtcur*xferrate(m)
         tmpa = g_soa - g_star(m)
         if (tmpa > 0.0_r8) then
            a_soa_tmp = a_soa(m) + beta(m)*tmpa
            sat(m) = g0_soa/(a_soa_tmp + a_opoa(m))
            g_star(m) = sat(m)*a_soa_tmp   ! this just needed for diagnostics
         end if
      end do
                                                                                                                                            
! step 2 - implicit in g_soa and semi-implicit in a_soa,
!    with g_star(m) calculated semi-implicitly
      tmpa = 0.0
      tmpb = 0.0
      do m = 1, ntot_soamode
         tmpa = tmpa + a_soa(m)/(1.0_r8 + beta(m)*sat(m))
         tmpb = tmpb + beta(m)/(1.0_r8 + beta(m)*sat(m))
      end do
                                                                                                                                            
      g_soa = (tot_soa - tmpa)/(1.0_r8 + tmpb)
      g_soa = max( 0.0_r8, g_soa )
      do m = 1, ntot_soamode
         a_soa(m) = (a_soa(m) + beta(m)*g_soa)/   &
                    (1.0_r8 + beta(m)*sat(m))
      end do
                                                                                                                                            
!     if (idiagss > 0) then
!        write(luna,'(i3,1p,20e10.2)') niter, tcur, dtcur, &
!           phi(:), g_star(:), a_soa(:), g_soa
!        write(luna,'(23x,1p,20e10.2)') &
!           sat(:), sat(:)*a_soa(:), a_opoa(:)
!     end if

!     if (niter > 9992000) then
!        write(luna,'(a)') '*** to many iterations'
!        exit
!     end if

      end do timeloop


      g_soa_tend = (g_soa - g_soa_in)/dtfull
      do m = 1, ntot_soamode
         a_soa_tend(m) = (a_soa(m) - a_soa_in(m))/dtfull
      end do


      return
      end subroutine modal_aero_soaexch

!----------------------------------------------------------------------


      subroutine modal_aero_gasaerexch_init

!-----------------------------------------------------------------------
!
! Purpose:
!    set do_adjust and do_aitken flags
!    create history fields for column tendencies associated with
!       modal_aero_calcsize
!
! Author: R. Easter
!
!-----------------------------------------------------------------------

use modal_aero_data
use modal_aero_rename

use abortutils, only   :    endrun
use cam_history, only  :   addfld, add_default, fieldname_len, phys_decomp
use constituents, only :  pcnst, cnst_get_ind, cnst_name
use spmd_utils, only   :    masterproc
use phys_control,only  : phys_getopts


implicit none

!-----------------------------------------------------------------------
! arguments

!-----------------------------------------------------------------------
! local
   integer  :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
   integer  :: jac
   integer  :: l, lsfrm, lstoo, lunout
   integer  :: l_so4g, l_nh4g, l_msag, l_soag
   integer  :: m, mfrm, mtoo
   integer  :: nsamefrm, nsametoo, nspec
   integer  :: n, nacc, nait

   logical  :: do_msag, do_nh4g, do_soag
   logical  :: dotend(pcnst), dotendqqcw(pcnst)

   real(r8) :: tmp1, tmp2

   character(len=fieldname_len)   :: tmpnamea, tmpnameb
   character(len=fieldname_len+3) :: fieldname
   character(128)                 :: long_name
   character(8)                   :: unit

   logical                        :: history_aerosol      ! Output the MAM aerosol tendencies

   !-----------------------------------------------------------------------
 
        call phys_getopts( history_aerosol_out        = history_aerosol   )
 
	lunout = 6
!
!   define "from mode" and "to mode" for primary carbon aging
!
!   skip (turn off) aging if either is absent, 
!      or if accum mode so4 is absent
!
	modefrm_pcage = -999888777
	modetoo_pcage = -999888777
	if ((modeptr_pcarbon <= 0) .or. (modeptr_accum <= 0)) goto 15000
	l = lptr_so4_a_amode(modeptr_accum)
	if ((l < 1) .or. (l > pcnst)) goto 15000

	modefrm_pcage = modeptr_pcarbon
	modetoo_pcage = modeptr_accum

!
!   define species involved in each primary carbon aging pairing
!	(include aerosol water)
!   
!
	mfrm = modefrm_pcage
	mtoo = modetoo_pcage

	nspec = 0
aa_iqfrm: do iqfrm = -1, nspec_amode(mfrm)

	    if (iqfrm == -1) then
		lsfrm = numptr_amode(mfrm)
		lstoo = numptr_amode(mtoo)
	    else if (iqfrm == 0) then
!   bypass transfer of aerosol water due to primary-carbon aging
		cycle aa_iqfrm
!               lsfrm = lwaterptr_amode(mfrm)
!               lstoo = lwaterptr_amode(mtoo)
	    else
		lsfrm = lmassptr_amode(iqfrm,mfrm)
		lstoo = 0
	    end if
	    if ((lsfrm < 1) .or. (lsfrm > pcnst)) cycle aa_iqfrm

	    if (lsfrm>0 .and. iqfrm>0 ) then
! find "too" species having same lspectype_amode as the "frm" species
		do iqtoo = 1, nspec_amode(mtoo)
		    if ( lspectype_amode(iqtoo,mtoo) .eq.   &
		         lspectype_amode(iqfrm,mfrm) ) then
			lstoo = lmassptr_amode(iqtoo,mtoo)
			exit
		    end if
		end do
	    end if

	    if ((lstoo < 1) .or. (lstoo > pcnst)) lstoo = 0
	    nspec = nspec + 1
	    lspecfrm_pcage(nspec) = lsfrm
	    lspectoo_pcage(nspec) = lstoo
	end do aa_iqfrm

	nspecfrm_pcage = nspec

!
!   output results
!
	if ( masterproc ) then

	write(lunout,9310)

	  mfrm = modefrm_pcage
	  mtoo = modetoo_pcage
	  write(lunout,9320) 1, mfrm, mtoo

	  do iq = 1, nspecfrm_pcage
	    lsfrm = lspecfrm_pcage(iq)
	    lstoo = lspectoo_pcage(iq)
	    if (lstoo .gt. 0) then
		write(lunout,9330) lsfrm, cnst_name(lsfrm),   &
      			lstoo, cnst_name(lstoo)
	    else
		write(lunout,9340) lsfrm, cnst_name(lsfrm)
	    end if
	  end do

	write(lunout,*)

	end if ! ( masterproc ) 

9310	format( / 'subr. modal_aero_gasaerexch_init - primary carbon aging pointers' )
9320	format( 'pair', i3, 5x, 'mode', i3, ' ---> mode', i3 )
9330	format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
9340	format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )


15000 continue

! set gas species indices
      call cnst_get_ind( 'H2SO4', l_so4g, .false. )
      call cnst_get_ind( 'NH3',   l_nh4g, .false. )
      call cnst_get_ind( 'MSA',   l_msag, .false. )
      call cnst_get_ind( 'SOAG',  l_soag, .false. )
      if ((l_so4g <= 0) .or. (l_so4g > pcnst)) then
         write( *, '(/a/a,2i7)' )   &
            '*** modal_aero_gasaerexch_init -- cannot find H2SO4 species',   &
            '    l_so4g=', l_so4g
         call endrun( 'modal_aero_gasaerexch_init error' )
      end if
      do_nh4g = .false.
      do_msag = .false.
      do_soag = .false.
      if ((l_nh4g > 0) .and. (l_nh4g <= pcnst)) do_nh4g = .true.
      if ((l_msag > 0) .and. (l_msag <= pcnst)) do_msag = .true.
      if ((l_soag > 0) .and. (l_soag <= pcnst)) do_soag = .true.

! set tendency flags
      dotend(:) = .false.
      dotend(l_so4g) = .true.
      if ( do_nh4g ) dotend(l_nh4g) = .true.
      if ( do_msag ) dotend(l_msag) = .true.
      if ( do_soag ) dotend(l_soag) = .true.
      do n = 1, ntot_amode
         l = lptr_so4_a_amode(n)
         if ((l > 0) .and. (l <= pcnst)) then
            dotend(l) = .true.
            if ( do_nh4g ) then
               l = lptr_nh4_a_amode(n)
               if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
            end if
         end if
         l = lptr_soa_a_amode(n)
         if ((l > 0) .and. (l <= pcnst)) then
            dotend(l) = .true.
         end if
      end do

      if (modefrm_pcage > 0) then
         do iq = 1, nspecfrm_pcage
            lsfrm = lspecfrm_pcage(iq)
            lstoo = lspectoo_pcage(iq)
            if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
               dotend(lsfrm) = .true.
               if ((lstoo > 0) .and. (lstoo <= pcnst)) then
                  dotend(lstoo) = .true.
               end if
            end if
         end do
      end if


!  define history fields for basic gas-aer exchange
!  and primary carbon aging from that
      do l = 1, pcnst
         if ( .not. dotend(l) ) cycle

         tmpnamea = cnst_name(l)
         fieldname = trim(tmpnamea) // '_sfgaex1'
         long_name = trim(tmpnamea) // ' gas-aerosol-exchange primary column tendency'
         unit = 'kg/m2/s'
         call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
         if ( history_aerosol ) then 
            call add_default( fieldname, 1, ' ' )
         endif
         if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit

      end do   ! l = ...


!  define history fields for aitken-->accum renaming
      dotend(:) = .false.
      dotendqqcw(:) = .false.
      do ipair = 1, npair_renamexf
         do iq = 1, nspecfrm_renamexf(ipair)
            lsfrm = lspecfrma_renamexf(iq,ipair)
            lstoo = lspectooa_renamexf(iq,ipair)
            if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
               dotend(lsfrm) = .true.
               if ((lstoo > 0) .and. (lstoo <= pcnst)) then
                  dotend(lstoo) = .true.
               end if
            end if

            lsfrm = lspecfrmc_renamexf(iq,ipair)
            lstoo = lspectooc_renamexf(iq,ipair)
            if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
               dotendqqcw(lsfrm) = .true.
               if ((lstoo > 0) .and. (lstoo <= pcnst)) then
                  dotendqqcw(lstoo) = .true.
               end if
            end if
         end do ! iq = ...
      end do ! ipair = ...

      do l = 1, pcnst
      do jac = 1, 2
         if (jac == 1) then
            if ( .not. dotend(l) ) cycle
            tmpnamea = cnst_name(l)
         else
            if ( .not. dotendqqcw(l) ) cycle
            tmpnamea = cnst_name_cw(l)
         end if

         fieldname = trim(tmpnamea) // '_sfgaex2'
         long_name = trim(tmpnamea) // ' gas-aerosol-exchange renaming column tendency'
         unit = 'kg/m2/s'
         if ((tmpnamea(1:3) == 'num') .or. &
             (tmpnamea(1:3) == 'NUM')) unit = '#/m2/s'
         call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
         if ( history_aerosol ) then 
            call add_default( fieldname, 1, ' ' )
         endif
         if ( masterproc ) write(*,'(3(a,3x))') 'gasaerexch addfld', fieldname, unit
      end do   ! jac = ...
      end do   ! l = ...


! calculate soa_equivso4_factor
! if do_soag == .false., then set it to zero as a safety measure
      soa_equivso4_factor = 0.0
      if ( do_soag ) then
         tmp1 = -1.0 ; tmp2 = -1.0
         do l = 1, ntot_aspectype
            if (specname_amode(l) == 's-organic') tmp1 = spechygro(l)
            if (specname_amode(l) == 'sulfate'  ) tmp2 = spechygro(l)
         end do
         if ((tmp1 > 0.0_r8) .and. (tmp2 > 0.0_r8)) then
            soa_equivso4_factor = tmp1/tmp2
         else
            write(*,'(a/a,1p,2e10.2)') &
               '*** subr modal_aero_gasaerexch_init', &
               '    cannot find hygros - tmp1/2 =', tmp1, tmp2
            call endrun()
         end if
      end if

      return
      end subroutine modal_aero_gasaerexch_init


!----------------------------------------------------------------------

end module modal_aero_gasaerexch


